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This paper describes a quantum Monte Carlo method capable of sampling the full density matrix 
of a many-particle system, thus granting access to arbitrary reduced density matrices and allowing 
expectation values of complicated non-local operators to be evaluated easily. The direct sampling 
of the density matrix also raises the possibility of calculating previously inaccessible entanglement 
measures. The algorithm closely resembles the recently introduced full configuration interaction 
quantum Monte Carlo method, but works all the way from infinite to zero temperature. We explain 
the theory underlying the method, describe the algorithm, and introduce an importance-sampling 
procedure to improve the stochastic efficiency. To demonstrate the potential of our approach, the 
energy and staggered magnetization of the isotropic antiferromagnetic Heisenberg model on small 
lattices and the concurrence of one-dimensional spin rings are compared to exact or well-established 
results. Finally, the nature of the sign problem in the method is investigated. 



I. INTRODUCTION 

Quantum Monte Carlo (QMC) methods are well es- 
tablished as vital tools in the study of complex many- 
body quantum systems, often providing highly accurate 
results. However, the use of QMC methods in the calcu- 
lation of quantum information measures has been some- 
what limited, largely due to their inability to sample re- 
duced density matrices involving more than a few parti- 
cles or lattice sites. 

Entanglement, well established as an important con- 
cept in quantum information theory, has recently become 
a subject of active research in the condensed matter com- 
munity, with various entanglement measures proving use- 
ful when studying how ground states change near quan- 
tum phase transitions In small systems the Lanczos 
method can be used to calculate exact entanglement mea- 
sures, and in one-dimensional systems the density matrix 
renormalization group method is applicable [H . However, 
the study of entanglement in large systems in more than 
one dimension is less straightforward. Relationships be- 
tween reduced density matrices and spin correlation func- 
tions [1] allow QMC methods to access some entangle- 
ment measures in certain situations Moreover, Sand- 
vik recently introduced a QMC method formulated in the 
valence-bond basis Q that allowed certain new entangle- 
ment entropy measures to be evaluated Hastings et 
al. showed that Renyi S2 entanglement entropy can also 
be calculated However, in general, the inability of 
QMC methods to directly access the density matrix and 
reduced density matrices of systems has hindered their 
use in this area. 

Projector QMC methods such as diffusion Monte Carlo 
[| (DMC) and Green's function Monte Carlo H, ^ 
(GFMC) grant access to zero-temperature properties by 
stochastically applying a projection operator to some ini- 
tial state such that a stochastic sampling of the ground- 



state wave function is eventually achieved. The fixed- 
node approximation often allows accurate results to 
be obtained in systems with sign problems, with the ac- 
curacy of such results depending on the quality of the 
nodal surface used. Finite-temperature methods take a 
different approach. Path-integral Monte Carlo (PIMC) 
methods express the partition function, Z = Tr(e~^^), 
as a sum over contributions from various paths in imag- 
inary time [ni . With an appropriate update procedure, 
the paths can be sampled with the correct probabilities, 
thus allowing finite-temperature expectation values to be 
obtained. The stochastic series expansion (SSE) method 
[T^ l has much in common with this approach but begins 
by performing a Taylor expansion of the density matrix. 
The resulting terms in the expansion are then sampled 
in a similar manner. These methods also allow access to 
ground-state properties. However, the sign problem is, 
in many cases, insurmountable at low temperatures. 

The full configuration interaction quantum Monte 
Carlo (FCIQMC) method recently introduced by Booth, 
Thom and Alavi [l3j is a projector method for studying 
zero-temperature properties, and, as such, has much in 
common with DMC and GFMC. However, unlike DMC 
and GFMC, where the sampling of the ground-state wave 
hmction is performed in real space, FCIQMC samples 
the wave function in a discrete basis. Crucially, no 
prior knowledge of the nodal structure of the ground- 
state wave function is required to reach the exact ground 
state. Rather, the sign problem manifests itself in the 
large but system-specific population of quantum Monte 
Carlo walkers required in order for the ground state of the 
Hamiltonian to emerge from the background noise. 
The system sizes accessible to FCIQMC are limited by 
the amount of memory available to store these walk- 
ers. However, the method has proven highly successful 
in many chemical systems, reducing the memory needed 
to achieve FCTquality results by several orders of mag- 
nitude This has led to much interest in this di- 
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rection and research into fundamental improvements and 
new applications of the algorithm continues [iB, [lO] • 

This article presents a closely-related QMC method, 
which we call density matrix quantum Monte Carlo 
(DMQMC). Like the path- integral and SSE methods, 
DMQMC allows finite-temperature results to be calcu- 
lated. However, it uses a projection approach to achieve 
this and thus has more in common with zero-temperature 
QMC methods. DMQMC was inspired by FCIQMC and 
shares many of its features. 

In DMQMC, rather than sampling the components of 
the wave function in a discrete basis, the elements of 
the density matrix are sampled instead. It is then a 
simple task to calculate expectation values of arbitrary 
operators, even those that do not commute with the 
Hamiltonian. Such expectation values arc usually dif- 
ficult to calculate using other QMC methods [2lj . More- 
over, the ability to directly sample the density matrix 
means that many quantum information measures are ac- 
cessible. These advantages cannot be expected to come 
without drawbacks, and, indeed, the systems to which we 
have successfully applied DMQMC to date are small by 
the standards of other finite-temperature methods. How- 
ever, the potential uses of a direct stochastic sampling of 
the density matrix are such that DMQMC deserves in- 
vestigation. This paper demonstrate the use of DMQMC 
by studying the isotropic antiferromagnetic Heisenberg 
model, but DMQMC is a general method and is applica- 
ble to both bosonic and fcrmionic systems. 

Section |ll] summarizes the FCIQMC method, setting 
the stage for the description of DMQMC in Section IIIII 
In Section IIVI an importance-sampling procedure is in- 
troduced. The DMQMC method is then applied to 
the isotropic antiferromagnetic Heisenberg model in Sec- 
tion |V] and used to calculate the energy and staggered 
magnetization of small square lattices and the concur- 
rence of onc-dimensional rings. The sign problem in 
DMQMC is investigated on the 4x4 triangular lattice 
in Section rvTl We discuss our results and offer some con- 
cluding remarks in Section [VIII Hartree atomic units are 
used throughout. 



II. FULL CONFIGURATION INTERACTION 
QUANTUM MONTE CARLO 

The DMQMC algorithm was formulated in analogy 
with the FCIQMC method. The two methods share 
many of the same features and it is often useful to com- 
pare them. We therefore begin by providing a brief sum- 
mary of the FCIQMC method; for more detailed discus- 
sions readers are referred to Refs. dH and dii). 

Consider the imaginary-time Schrodingcr equation 



dr 



The general solution to this equation is 



|*(r)) 
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(1) 



(2) 



for some initial wave function |'0(r = 0)). If the initial 
wave function has a non-zero ground-state component, 
co(0), it is easy to see that |^'(t)) will become propor- 
tional to the ground state in the limit of large r. 



|«'(t^c»)) =co(0)e-^^«|£^o> 



(3) 



where \Eo) is the ground state and Eq the ground-state 
energy. The factor of e~^°'^ can be removed by choosing 
the zero of energy such that = Q. In practice, since 
i?o is usually unknown, we solve 



dr 



= -{H - si) I*) = f I*) 



(4) 



where we have defined T — —{H — Si). The energy shift 
S is adjusted slowly during the simulation to keep the 
normalization approximately constant. The long-time 
average of S gives another measure of the ground-state 
energy. 

The above theory is common to all projector meth- 
ods; the difference is how they achieve the evolution to 
the ground state. FCIQMC works in a discrete basis 
of kcts \Xi)^ which are normally Slater determinants for 
fcrmions or permanents for bosons. The components of 
the wave function in this basis are represented stochas- 
tically by a collection of markers. Following Anderson 
[2^ [23I , we refer to these markers as "psips" . Each psip 
has an associated sign, g = ±1, which we refer to as its 
"charge", and resides on a particular basis state \Xi) (or 
on "site i"). The net psip charge on a basis state is in- 
terpreted as the amplitude of that state in the expansion 
of the wave function. At any point in a simulation, the 
distribution of psip charges does not need to provide an 
accurate representation of the wave function. Rather, 
the FCIQMC method only requires that the expectation 
values of the site charges represent the ground state . 
Thus, at any point in the simulation, the memory re- 
quired to sample the wave function may be many orders 
of magnitude smaller than the memory required to store 
the whole state. 

Booth and co-workers [l^ introduced an algorithm 
to evolve the population of psips according to the 
imaginary-time Schrodingcr equation. This can be sum- 
marized as follows. For each time step At we loop over 
the entire population of psips and perform the following 
steps: 

1. Spawning: Allow a psip with charge qi on site i to 
attempt to spawn onto a connected site j, where 
Tij 7^ and i 7^ 7, with probability |Tji|Ar. If the 
spawning attempt is successful, a psip is born at 
site i with charge qj = iiign[Tji)qi. 

2. Diagonal death/ cloning: Each psip has the chance 
to either clone or die with probability |Tjj|AT. The 
consequence of a successful death/cloning event de- 
pends on the sign of the diagonal matrix element: 
if > the psip is cloned; otherwise the psip is 
removed from the simulation. 
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3. Annihilation: Pairs of psips on the same site with 
opposite charges cancel out ("annihilate") and are 
removed from the simulation, leaving a population 
of only a single charge type on each site. 

The FCIQMC algorithm samples the solution of a first- 
order Euler finite-difference approximation to Eq. ([4]). 
Hence, the distribution of psips gives a stochastic repre- 
sentation of the wave function and, as r — )■ (X), the psips 
settle down to sample the ground-state wave function. 

The psip annihilation step does not alter the total 
charge on a site. However, it has been shown to be 
vital in order for the true ground-state wave function 
to emerge in systems with sign problems [13 ]. Similar 
and more complex walker cancellation mechanisms have 
been attempted in continuum DMC and GFMC calcu- 
lations with less success p3 - l29j . Walker cancellation 
in FCIQMC works better because psips are more likely 
to encounter each other in a discrete and finite Hilbert 
space. 

The ability of FCIQMC to tackle the sign problem 
through the annihilation step is perhaps its most signifi- 
cant advantage. Annihilation leads to the characteristic 
population dynamics and allows the true ground state 
to emerge without any knowledge of the nodal structure 
of the ground-state wave function [ij]. At the start of 
an FCIQMC simulation, the shift is held constant at a 
value large enough to ensure that the psip population 
grows exponentially. Eventually, when a system-specific 
population of psips is reached, the annihilation rate be- 
comes equal to the spawning rate and the population 
spontaneously plateaus. During this plateau period the 
ground state of the correct Hamiltonian emerges, after 
which the population begins to grow again. The popula- 
tion can then be controlled by varying the shift. The psip 
population at the plateau must be a small fraction of the 
number of basis states in the Hilbert space in order for 
FCIQMC to be more (memory) efficient than an exact 
diagonalization. 

The ground-state energy in FCIQMC can be calculated 
using a "projected estimator" based on the the expression 



where = {X,\0\Xj). Although 



Xi, the ex- 



Eo 



T.i.j^lHtjX] 



(5) 



Here |*^) is an appropriate trial function with compo- 
nents ijjj in the many-particle basis, jA^i), and Xi is a 
component of the exact ground state in this basis. The 
ground-state energy is obtained by averaging the nu- 
merator and denominator separately, with the total psip 
charge on each site, q\°^ , used in the place of correspond- 
ing exact amplitude, Xi- 

Calculating the ground-state expectation value of an 
operator O that does not commute with the Hamiltonian 
is more difficult because a projected estimator cannot be 
used. Instead, assuming that \E^) is real, it is necessary 
to evaluate 



(0) = 



(go|Q|go 
{Eo\Eo) 



(6) 



pectation value of a product is not equal to the product 
of the expectation values: ^ XiXj- This means 

that XiXj cannot be obtained by averaging the products 
of the instantaneous psip weights. One could in princi- 
ple average 9'°* over many iterations to obtain Xi before 
multiplying Xi ^-nd Xjj but this would involve storing a 
number for every basis function, which is unfeasible due 
to memory limitations. 

This problem is not easy to overcome and there is cur- 
rently no way to evaluate general expectation values ex- 
actly within the FCIQMC framework. Indeed, the calcu- 
lation of general expectation values in other Monte Carlo 
methods is often a difficult task. 



III. DENSITY MATRIX QUANTUM MONTE 
CARLO 



We now show how an FCIQMC-like dynamics can be 
used to sample both finite-temperature and ground-state 
density matrices. We first consider the thermal density 
matrix and how it can be evolved as a function of inverse 
temperature by solving the symmetrized Bloch equation. 
We then draw upon analogies with FCIQMC to formulate 
the DMQMC algorithm before discussing the calculation 
of estimators for a general quantum mechanical observ- 
able. This section ends with an explanation of how to 
sample a reduced density matrix in order to calculate 
estimators of entanglement measures. 



A. Theory 

Since the psip population (and hence the normaliza- 
tion) varies during a quantum Monte Carlo simulation, 
it is convenient to work with the unnormalized thermal 
density matrix 



(7) 



where H is the Hamiltonian operator and /3 = l/ksT is 
the inverse temperature. The canonical partition func- 
tion Z{l3) is given by: 



Z(/3) = Tr(p(/?)). 



(8) 



Differentiating p(j3) with respect to (3 shows that it obeys 
both the Bloch equation. 



dp /y. 

dp = 

and the symmetrized Bloch equation. 



(9) 



(10) 



The symmetrized version turns out to be more use- 
ful for our purposes. Consider the general solution to 
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Eq. (nnD, 

p{P)=e-i"p{P^i))e-i", (11) 

and let \Ei) denote an energy eigenstate with energy Ei. 
Without loss of generality, we assume that i?o = 0. By 
expanding p{P = 0) as 

p{P^Q)=Y,p,Mm{E,\ (12) 

and applying Eq. pT|) it is seen that p(/3) becomes pro- 
portional to the ground-state density matrix as /? — >■ oo: 

/5(/3^(») =poo(0)|So)(£^o|. (13) 

A similar analysis shows that the ground-state density 
matrix is not reached if a general initial density ma- 
trix is propagated using the unsymmetrized Bloch equa- 
tion. Hence, in cases where ground-state properties are 
desired, Eq. (fTO|) is more useful. If we wish to sample 
thermal properties at any finite temperature T, we must 
also impose the correct initial condition on Eq. (fTO)) . The 
most convenient is at /3 = 0, where 

p(/3 = 0) = i. (14) 

Eqs. (|10[) and p4p provide everything required to 
stochastically sample the thermal density matrix. As 
in FCIQMC we introduce a collection of psips, which in 
DMQMC occupy not basis states \Xi) but basis operators 
\Xi) {Xj\ of the expansion of the density operator /3(/?) 
in some convenient basis set. At the start of a simulation 
the psips are distributed randomly along the diagonal of 
the infinite-temperature density matrix Pij{P = 0) = (5^. 
The psips then evolve under a set of rules so that they 
sample the matrix elements Pij{l3) in the chosen basis 
at each iteration. In order to prevent the population of 
psips from exploding we introduce an energy shift, 5, and 
let H ^ H - si as in FCIQMC. Eq. ^ thus becomes 

|4(TP + Pf), (15) 
where T = —{H — Si) is the "update matrix". 

B. Algorithm 

We use the following algorithm to evolve a collection of 
psips according to Eq. (fT5|) . For a single step in inverse 
temperature of A/3, we loop over the entire population of 
psips and perform the following steps: 

1. Spawning along columns of the density matrix: Al- 
low a psip with charge Qij on site to attempt 
to spawn onto connected sites {k,j), where Tik ^ 
and i k, with probability i|Ti;;|A/3. If the spawn- 
ing attempt is successful, a psip is born at {k,j) 
with charge qtj = sigii(Tik)qij. 



2. Spawning along rows of the density matrix: Simi- 
larly, allow a psip with charge qij on site {i,j) to 
attempt to spawn onto connected sites (i, k) with 
probability i|Tjfc|A^. If the spawning attempt is 
successful, a psip is born at site {i, k) with charge 
Qik = sign{Tjk)qtj. 

3. Diagonal death/cloning: If there is a psip on site 
(i, j) and Ta +Tjj < 0, then with probability pd = 
'^\Tii + Tjj\Af3 that psip is killed and removed from 
the simulation; if Tu + Tjj > 0, the psip is cloned 
(i.e., a new psip is created on the same site and 
with the same charge) with probability pd- 

4. Annihilation: Pairs of psips inhabiting the same 
site (i,j) with opposite charges annihilate and are 
removed from the simulation, leaving a population 
of only a single charge type on each site. 

The first three steps describe a stochastic algorithm to 
sample the solution of a first-order Euler finite-difference 
approximation to Eq. pO[) . The distribution of psip 
charges a.t /3 + A/3 is thus proportional to the density 
matrix at this inverse temperature, provided that the dis- 
tribution of charges at /3 was correct. 

The DMQMC method shares many similarities with 
FCIQMC, namely: 

• Annihilation does not alter the expected (normal- 
ized) psip distribution but serves to overcome the 
sign problem 

• The underlying finite-difference approximation is 
stable if < A/3 < 2/(£',„ax — Eq), where -Emax 
is the largest eigenvalue of the Hamiltonian matrix 
[l^ . This is a sufficient condition to ensure cor- 
rect projection onto the exact ground state, but 
the finite value of A/3 leads to an error of C'(A/3) 
in the density matrix at temperatures greater than 
zero. It is therefore necessary to check that finite- 
temperature results have converged with respect to 
A/3. 

• The familiar shift-update algorithm already used 
in DMC HO] and FCIQMC [ill simulations is em- 
ployed to modify S and thus to control the popu- 
lation. The shift, S, is adjusted according to 

5(/3+AA/3) = SiP) 

C . / jVp(/3 + AA/3) \ (16) 
AA/3 Np{l3) J ' 

where A is the number of /3-steps between shift 
updates, C is a shift damping parameter, and 
7Vp(/3) is the total number of psips at the inverse- 
temperature (3. During simulations, C, is chosen 
carefully to prevent large fluctuations in S. 

• Rather than attempting all possible spawning 
events from a psip on a given site, it is compu- 
tationally efficient to attempt just one (or a small 
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number) of the many possible spawning events and 
reweight the acceptance probabihties accordingly 

m- 

• The algorithm is highly parallelizable — only the 
annihilation step requires communication between 
CPU cores and so using a large psip population is 
a viable option. 

As with all projection methods, the ground state is 
approached as /3 oo. Once convergence to the ground 
state has been attained the density matrix no longer de- 
pends on (3 (remember that S is chosen to keep the nor- 
malization fixed) and the statistical errors of measured 
ground-state expectation values can easily be reduced by 
averaging over imaginary time. A single simulation is 
therefore sufficient to obtain accurate ground-state ex- 
pectation values. The estimation of finite-temperature 
properties is more difficult because the inverse tempera- 
ture /? changes continuously as the simulation progresses; 
it is not possible to hold the simulation at a specific tem- 
perature whilst statistics are accumulated. 

With access to a stochastic sampling of the unnormal- 
ized density matrix at a given temperature, the expecta- 
tion value of any quantum mechanical observable, O, at 
that temperature can be calculated using 



(O) 



Tr(/5) 



77*°' (9 



(17) 



The numerator and denominator must be sampled and 
averaged separately at each temperature to achieve the 
desired statistically accuracy. Hence a simulation in- 
volves repeatedly evolving the density matrix from /3 = 
to some chosen maximum value of /?, a process we call a 
"/3-loop" . This also gives us the freedom to allow Eq. p4|) 
to be satisfied exactly only on average, which greatly re- 
duces the memory demands. Each /3-loop is initialized 
with a different random number seed and with psips ran- 
domly distributed with uniform probability along the di- 
agonal of the density matrix. Statistics are accumulated 
over a sufficiently large number of /3-loops in order to 
obtain the desired statistical accuracy. As each /3-loop 
is completely independent, performing multiple /3-loops 
is an embarrassingly parallel task and each /3-loop gives 
statistically independent data points. 

A primary concern is that averages over an impracti- 
cally large number of /3-loops may be required. Indeed, 
it is not uncommon to have to average over 0(10^) itera- 
tions in ground-state calculations. In practice, however, 
reaching a sufficiently large number of /3-loops does not 
appear to be a problem. The crucial difference from the 
ground-state case is that each /3-loop provides statisti- 
cally independent data whereas ground-state calculations 
need to overcome inherent correlations between nearby 
iterations. Furthermore, as will be shown, the quality 
of estimates obtained using a given number of /3-loops is 
better at high temperature than at low temperature. An- 
other advantage is that a single simulation can provide 
data across the entire temperature range being studied. 



The energy shift, S, is varied in order to control the 
normalization. This, however, introduces a bias: if the 
shift is varied by AS, then the psip population is effec- 
tively altered by a factor e'^^^'^. As a result, configura- 
tions in which the average energy of the psip population 
happens to be less negative than usual are effectively 
given too large a weight. This problem is particularly 
severe in DMQMC because the results at each tempera- 
ture are obtained by averaging over separate /3-loops, and 
thus, for each quantity contributing to a given estimator, 
the corresponding shift (and hence population) profiles 
can be very different. 

The population bias can be greatly reduced using a 
method suggested by Umrigar et al. in the context of 
DMC (30| . in which each sampled quantity proportional 
to the psip population is multiplied by the factor 



B~l 



n(/3,i3)= n 



-S(,3-mA/3) 



(18) 



where B is some chosen number of factors and B = 
mm{-^,B). By multiplying by this factor, we remove 

the last B factors of e^^^"^ introduced by varying the 
shift. As S — > oo the population control bias should 
be completely removed. The population control bias can 
also be reduced by using a larger population of psips. 

Another concern is that the number of elements in 
the density matrix is the square of the dimension of the 
Hilbert space of many-particle states. At first glance it 
seems that, without dramatically increasing the number 
of psips, the density matrix would be very poorly sam- 
pled. However, the dimension of the Hilbert space, D, 
rises exponentially with the number of particles or sites, 
N, so that D (X e"^ and thus cx e"^^^'). The dou- 
bling of the exponent implies that a DMQMC simula- 
tion for an Af/2-site lattice model requires approximately 
the same number of psips as an TV-site FCIQMC simu- 
lation to achieve the same sampling quality. Moreover, 
DMQMC estimators for operators that do not commute 
with the Hamiltonian often have significantly smaller 
variance than the forward-walking or other estimators 
required to evaluate the ground-state expectation values 
of such operators in FCIQMC. 



C. Entanglement measures 

It straightforward to obtain a stochastic representation 
of any reduced density matrix (RDM) from a stochastic 
representation of the full density matrix. Reduced den- 
sity matrices are required to calculate entanglement mea- 
sures such as the concurrence described in Appendix [X] 

An RDM can be sampled by "tracing out" unwanted 
psips. Consider a composite quantum system C, which 
can be partitioned into two subsystems A and B so that 
He = Ha (^Hb, where Ha (He) denotes the Hilbert 
space of system A (B). The RDM pA that describes 
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sublattice A is defined by taking the partial trace of the 
fuU density matrix, pc, over ah the sites on sublattice B: 

PA^TvBipc)- (19) 

Our implementation of DMQMC represents the many- 
particle basis functions as bit strings [sif . where each 
bit refers to the state of a single spin. To evaluate pA 
we construct a mask. Ib, which only has bits set which 
correspond to spins in subsystem B. The density 
matrix element of pc contributes to pA if the result of 
the logical AND operation of the i string with the Ib mask 
is identical to that of the j string. The corresponding 
element in the RDM can be found by taking AND with an 
analogous I a mask. 

We calculate the desired entanglement measure from 
the unnormalized stochastic RDM along with its trace at 
the desired temperature, as discussed in Appendix 1X1 for 
the case of concurrence. The entanglement measure and 
trace are averaged over multiple /3-loops, or, for ground- 
state properties, over many iterations starting from a suf- 
ficiently low temperature. 

Depending on the specific Hamiltonian, psip popula- 
tion and temperature, it may be that only a few psips 
survive the masking process to contribute to the RDM. 
In such cases it is necessary to increase the psip popula- 
tion or accumulate the RDM over several /3-loops before 
making a single estimate. The use of translational and 
other symmetries can improve the quality of the sampling 
of an RDM, although symmetry was not exploited in this 
work. 

Any reduced density matrix can be obtained using 
DMQMC, but this does not necessarily allow any entan- 
glement measure to be calculated. For example, calcu- 
lating the von Neumann entropy involves taking the log- 
arithm of the RDM, which requires the full RDM to be 
represented and diagonalized. This places a limit on the 
size of the reduced subsystem but not on the full system. 
Fortunately, many useful entanglement measures involve 
subsystems containing only a small number of spins. 

IV. IMPORTANCE SAMPLING 

This article considers the application of DMQMC to 
the 5 = 1/2 antiferromagnetic Heiscnberg model, 

H^JY,S^■S,, (20) 

where J > and the (i, j) implies that the summation is 
over nearest-neighbor pairs of spins only. Periodic bound- 
ary conditions are applied. We work in the standard basis 
set where the many-spin states are tensor products of the 
one-spin eigenstates, |t) and ||) , of the Sz operator. 

The algorithm described in section II allows a stochas- 
tic sampling of the exact finite-temperature density ma- 
trix, assuming A/3 is sufficiently small. However, for 
the antiferromagnetic Heiscnberg model, the sampling 



method described is found to be insufficient for all but 
the smallest systems. The ground-state wave function is 
highly delocalized over the states in the basis set, and 
thus the ground-state density matrix, p = \Eq) {Eq\, has 
many off-diagonal elements with magnitudes comparable 
to the diagonal elements. As a result, when the den- 
sity matrix is sampled via the DMQMC algorithm, only 
a small fraction of psips reside on or near diagonal ele- 
ments. Estimators for the expectation values of most op- 
erators of interest only receive contributions from psips 
on or near the diagonal, so very few psips contribute and 
these estimators suffer from large statistical errors at low 
temperatures. Importance sampling can greatly improve 
the sampling quality. 

We start by defining the excitation level between basis 
states \Xi) and \Xj) to be the smallest number of pairs 
of opposite spins that must be flipped in order to reach 
\Xi) from \Xj). For the Heiscnberg model, Eq. (PP)) . the 
excitation level can change by at most ±1 in a single 
application of the Hamiltonian. 

One straightforward way to improve the quality of 
sampling is to reduce the probability of psips spawn- 
ing far from the diagonal of the density matrix. Psips 
that do reside on higher excitation levels are given a cor- 
respondingly larger weight, so that expectation values 
of operators are unchanged. However, the increase in 
the population of low- weight psips near the diagonal of 
the density matrix reduces the stochastic error in near- 
diagonal expectation values. 

With this motivation in mind, we define the following 
importance-sampling procedure; a more rigorous formu- 
lation is given in Appendix |B] Every time a psip on exci- 
tation level 7 attempts to spawn a new psip on excitation 
level (5, the probability of successful spawning is altered 
by a factor P^g. Thus, if a psip on a diagonal element 
attempts to spawn a new psip onto the first excitation 
level, the probability of successful spawning is altered by 
a factor Pqi, where Poi < 1- The first excitation level 
will thus be occupied by Pqi times as many psips as it 
would have been with unaltered spawning probabilities. 
The reduced (relative) population must be accounted for 
by giving all psips in the first excitation level a weight of 
Wi ~ 1/^01 when evaluating estimators. The number of 
spawning events from the first to the second excitation 
level is also altered by a factor Pqi (due to the reduced 
number of psips at level 1) and the chance of success- 
ful spawning for each attempt is multiplied by the factor 
Pi2. Hence, the weight given to the second excitation 
level must be W2 = ^/PqiPi2- Furthermore, since there 
are Pqi times as many psips on the first excitation level, 
the probability of spawning from the first excitation level 
to the diagonal elements must be enhanced by a factor 
1/Poi in order to achieve consistent spawning dynamics. 

In general P-^g — 1 / Pg-y and the weight given to a psip 
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FIG. 1. (Color online.) The imaginary-time dependence of 
the fractions of psips on the first four excitation levels dur- 
ing a DMQMC simulation of a 4x4 square antiferromagnetic 
Ifeisenberg lattice 32] . A single /3-loop was carried out using 
an initial population of 10^ psips on the diagonal elements. 
The results shown in (a) were obtained without importance 
sampling. The results shown in (b) were obtained using im- 
portance sampling to ensure that every excitation level had a 
roughly equal number of psips at zero temperature. 



on excitation level 7 is 



(5—7 



The estimator for an expectation value, (O), is thus 



(21) 



(0> = 



J2^ h 



(22) 



where E{i,j) is the excitation level between \Xi) and 
\Xj), and pij is the importance-sampled density ma- 
trix. The expected value of the importance-sampled psip 
charge on site (i, j) is proportional to pij. 

There is clearly some freedom in the numerical val- 
ues chosen for the factors P^g. In this study we adjust 
P^s such that all excitation levels have similar psip pop- 
ulations in the ground state. Whilst this choice is not 
necessarily optimal, it successfully increases the quality 
of ground-state estimates by many orders of magnitude 
whilst still allowing the entire density matrix to be sam- 
pled. Figures [IK a) and[T](b) show the fractions of psips on 
the first four excitation levels for the 4x4 square Heisen- 
berg lattice as obtained without and with importance 
sampling. Both simulations used an initial population 
of 10^ psips and a single /3-loop. It is clear that impor- 
tance sampling greatly assists in keeping a non-negligible 
population on the lower excitation levels. 

In practice, the values of P^s are very small for small 
7 and S and decrease in magnitude as the lattice size 



increases. In order to avoid an unnecessarily large sup- 
pression of psip spawning from the diagonal at high tem- 
peratures, we introduce the weights and probabilities 
P-yS gradually from (3 = until they reach their desired 
final values at /3targct- For (3 < /3target, the weight at ex- 
citation level 7 is set equal to W^^^""''°'; for /3 > /Stargct, 
the weight is held constant. The value of /Stargot is cho- 
sen to be the imaginary time at which the ground state is 
deemed to have been reached. However, the value of this 
parameter is not critical for the quality of the sampling. 



V. RESULTS 

For this first study of DMQMC we focus primarily on 
reproducing some exact and well-established results. It 
is generally the case that QMC simulations of the an- 
tiferromagnetic Heisenberg model on square lattices of 
even dimensions do not suffer from the sign problem. In 
the FCIQMC and DMQMC methods this is because all 
psips spawned on any specific site have the same sign 
and no annihilation occurs (T^ . This system is therefore 
particularly convenient to study. 

In order to calculate real finite-temperature properties, 
it is necessary to include contributions from all Ms (total 
spin) subspaces. The Heisenberg Hamiltonian conserves 
A/5, so different subspaces can be studied using separate 
simulations, which is an embarrassingly parallel compu- 
tational task. Combining results for different values of 
Ms is straightforward but requires additional calcula- 
tions and reveals nothing of interest about the perfor- 
mance of DMQMC. We therefore present results for the 
Ms = subspace only. 

We first consider the example of the 4x4 lattice, which 
is small enough that an exact diagonalization can be per- 
formed, thereby allowing a direct check of our DMQMC 
results. Figure ^a.) shows the energy as a function of 
temperature using an initial population of 100 psips at 
the start of each /3-loop and accumulating statistics over 
1000 such loops. The shift was allowed to vary through- 
out the simulation and there were typically 700-800 psips 
in the simulation by the end of each /3-loop. No impor- 
tance sampling was applied, and so statistical flucuations 
increase with inverse temperature as explained in Sec. lIVI 
The agreement with the exact results is very good. In- 
creasing the initial population to 10^ psips (Fig.[2Kb)) re- 
duces the statistical errors such that the agreement with 
the exact results is essentially perfect. 

The square of the staggered magnetization is repre- 
sented by the operator 

= M • M, with M = ^(-l)^-+''-S',, (23) 

i 

where Xi and Xj denote the coordinates of the square lat- 
tice. This operator does not commute with the Hamilto- 
nian. Its expectation value, (M^), is plotted as a func- 
tion of P in Fig. [3] for an 8x8 lattice; a ground-state 
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FIG. 2. (Color online.) The exact and DMQMC ener- 
gies for the 4x4 antiferromagnetic Heisenberg lattice (Hilbert 
space dimension D ~ 1.29 x 10^) with periodic boundary 
conditions Hi]- The /? step size, A/3, obeys JiVA^ = 0.1, 
where A'^ = 16, and results are shown every 30 iterations. 
Each simulation consisted of 1000 /3-loops. For the simula- 
tion reported in (a), 100 psips were introduced at the start 
of each /3-loop; for (b), 10'' psips were introduced. The error 
bars in (b) are smaller than the size of the markers. 
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FIG. 3. (Color online.) The square of the staggered mag- 
netization for an 8x8 antiferromagnetic Heisenberg lattice 
{D « 1.83 X 10^**) using 1.4 x lO'' psips and 143 ^-loopsjli]. 
The importance-sampling procedure described in the main 
text was applied. The DMQMC value of (AP) is plotted ev- 
ery 50th iteration. Error bars, where not visible, are smaller 
than the size of the marker. The grou nd- state value of (M^) 
obtained from an SSE simulation '33\ is plotted for compari- 



value obtained using the SSE method [33| is also shown. 
The dimension of the Ms = Hilbert space of an 8x8 
Heisenberg Hamiltonian is approximately 1.83 x 10^®, so 
the density matrix sampled in this simulation has ap- 
proximately 3.36 X lO'^^ elements. 

We also considered the application of DMQMC to 
ground-state calculations, where statistics are accumu- 
lated over a single /3-loop after the ground state has 



been reached. A relatively modest simulation of a 6x6 
lattice using 6.5 x 10^ psips gave a ground-state en- 
ergy of -0.67888(24)JiV and a value of (Af^) equal to 
0.20985(7), where errors were obtained via a bloc king 
analysis [3J|. These values agree with exact results [33 
of -0.678872 J7V and 0.20983 for the ground-state energy 
and staggered magnetization, respectively, to within sta- 
tistical errors [36|. Remarkably, despite the fact that 
[AP,!!] 7^ 0, the statistical error in the value of (M^) 
was smaller than that of the energy. It was not necessary 
to reweight these results using Eq. (|18|) as the psip popu- 
lation was large enough to render such biases negligible. 

The ground-state concurrence, Cgs, for neighboring 
spins (qubits) on an antiferromagnetic Heisenberg ring 
was studied by Woottcrs and O'Connor in 2001 [37| . 
They showed that, for an even number of spins, Cgs has 
a simple relationship with the ground-state energy. 



Ca 



1) 



(24) 



The exact results calculated from this formula provide a 
useful test of our DMQMC estimates of Cgs, which are 
obtained from the sampled reduced density matrix (see 
Appendix |^ . 

The DMQMC estimates of Cgs are presented in Table[Tl 
along with Woottcrs and O'Connor's exact analytic val- 
ues for up to = 10 sites. A ring with A^ = 36 sites was 
also studied, using FCIQMC to calculate Eq and then 
Eq. p4|) to obtain a value of Cgs for comparison with the 
DMQMC results. For lengths up to A^ = 10, the calcu- 
lation of the concurrence can easily be carried out using 
other methods and provides a straightforward test of the 
DMQMC algorithm. The A^ = 36 chain is far from triv- 
ial, and it is promising that such accurate results can be 
obtained. 

We note that the DMQMC estimates of Cgs were ob- 
tained by sampling the reduced density matrix for a sin- 
gle pair of spins. Due to translational invariance, these 
estimates could have been improved by sampling the re- 
duced density matrix for every neighboring pair and com- 
bining the results. 



VI. THE SIGN PROBLEM IN DMQMC 

As discussed in Sec. [U the annihilation step in 
FCIQMC leads to the characteristic population dynam- 
ics, whereby the sign problem can only be overcome once 
a critical psip population (the 'plateau') has been ex- 
ceeded. Since the annihilation steps in FCIQMC and 
DMQMC are identical, we would expect DMQMC to 
possess similar population dynamics. The annihilation 
rate for a given psip population in DMQMC will be sig- 
nificantly smaller than in FCIQMC because the number 
of density matrix elements is the square of the number 
of basis functions. As such, it is to be expected that the 
plateau height for a given Hamiltonian will be higher in 
DMQMC than in FCIQMC. 
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AT 


Eo/N 


Exact [37] DMQMC Np 




n 


4 


-0.5000 


0.5000 


0.5005(4) 2.5x103 


250 


6 


6 


-0.4671 


0.4343 


0.4342(5) 2.5x10^ 


1x10^ 


20 


8 


-0.4564 


0.4128 


0.4129(5) 1x10** 


1x10^ 


70 


10 


-0.4515 


0.4031 


0.4031(4) 4x10^ 


1x10^ 


252 


36 


-0.44374(2) 0.38748(4) 0.3873(8) IxlO'^ 


12 


9.08 X 10^ 



TABLE I. DMQMC estimates of the ground-state concur- 
rence, Cgs , for antiferromagnetic spin rings containing TV sites 
in the absence of an external magnetic field. Each ring cor- 
reponds to a Hilbert space of D basis functions. The density 
matrix was sampled using Np psips, and statistics accumu- 
lated over Ni /3-loops. Exact results for the energy and con- 
currence for rings up to A'^ = 10 are taken from Ref. (|37l). 
An FCIQMC calculation was used to find the energy of the 



mined using Eq. (|24|) . The full density matrix of the A'" = 36 
chain has approximately = 8.24 x 10^^ elements. 



To investigate this, we considered the antiferromag- 
netic Heisenberg model on a 4x4 triangular lattice, for 
which an exact diagonalization is easily carried out. The 
triangular lattice is the archetypal example of a frus- 
trated lattice and has a severe sign problem. We carried 
out a single /?-loop and allowed the population to grow 
with a fixed shift whilst simultaneously investigating the 
accuracy of the DMQMC energy estimate. Figure [H^a) 
demonstrates that the population plateaus as expected. 
Figure |4ljb) shows the accuracy of the energy estimate 
throughout the plateau period. 

At high temperatures accurate results are obtained. 
This is also the case in other finite-temperature methods, 
where it is generally found that the sign problem is less 
severe for small /3; indeed, there is no sign problem at all 
at infinite temperature. However, at lower temperatures 
the energy estimates suffer large fluctuations, and it is 
not until the population exceeds the plateau that a good 
agreement with exact results is once again obtained. The 
plateau height occurs at ~ 1.35 x 10^ psips. For compari- 
son, the plateau population in an FCIQMC calculation of 
the same system is at ~ 1.25 x 10^ psips. We find (in this 
case) that the DMQMC plateau height is approximately 
the square of the plateau height in FCIQMC, matching 
the increase in the size of the space being sampled. It 
is possible to obtain accurate results for the entire tem- 
perature range simply by starting the simulation with 
an initial population greater than that of the plateau. 
Moreover, even if the plateau cannot be reached due to 
memory restrictions, one can nevertheless systematically 
reach lower temperatures by increasing the population of 
psips. 

The sign problem in the antiferromagnetic Heisenberg 
model on a triangular lattice is severe in both FCIQMC 



2; 




and DMQMC [IJ]. However, the efficiency of the anni- 
hilation procedure in FCIQMC varies substant ially with 
the system studied and with the basis set used [ij. Fur- 



FIG. 4. (Color online.) Population dynamics and energy es- 
timate for the antiferromagnetic Heisenberg model on a 4x4 
triangular lattice with periodic boundary conditions. A fixed 
shift was used and a single /3-loop performed. Panel (a) shows 
the emergence of a plateau in the psip population and (b) 
shows the exact energy and the DMQMC estimate of the en- 
ergy as functions of p. Once the plateau has been exited, 
the DMQMC energy is in good agreement with the exact re- 
sults. The severity of the sign problem increases with inverse 
temperature. 



thermore, the initiator app roximation (i-FClQMC) pro- 
posed by Cleland et al. [la| , has been shown to reduce the 
memory requirements for FCIQMC calculations by sev- 
eral orders of magnitude in many cases. In i-FClQMC, 
spawning events onto previously unoccupied sites are for- 
bidden unless the psip population of the parent site ex- 
ceeds a threshold. This increases the annihilation rate 
relative to FCIQMC and ameliorates the sign problem at 
the expense of introducing a systematically improvable 
approximation. We have yet to investigate the DMQMC 
equivalent of the initiator method. It is expected that 
ground-state properties will be available, but it is not 
yet clear to what extent the modified spawning will af- 
fect finite-temperature results. 



VII. DISCUSSION 

This article has described DMQMC, a quantum Monte 
Carlo method that allows direct sampling of the finite- 
temperature and ground-state density matrices in a dis- 
crete basis. The validity of the method has been veri- 
fied by reproducing exact and well-established results for 
some small systems, including the calculation of the con- 
currence of one-dimensional spin rings by directly sam- 
pling reduced density matrices. In all cases investigated, 
DMQMC has proved capable and accurate. 

The introduction of an importance-sampling procedure 
allows larger lattices to be investigated — the largest 
system that we have simulated successfully to date is a 
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10x10 antiferromagnetic Hcisenberg model. Larger sys- 
tems could easily be tackled at high temperatures, but it 
is unlikely that our very simple approach to importance 
sampling will allow simulations of lattices of more than 
lOx 10 sites over the entire range of temperatures from in- 
finity to zero. If a better importance-sampling procedure 
can be devised, as we believe likely, many more avenues 
will be opened. 

Like FCIQMC, DMQMC uses an annihilation pro- 
cedure that allows it to overcome the sign problem if 
a system-specific population of psips is reached. In 
FCIQMC the sign problem can often be ameliorated, 
for example by changing to a more appropriate basis set 
[3, 113 or by applying the initiator approximation [l^ . 
Given the similarities between DMQMC and FCIQMC, 
it is likely that these ideas and future developments in 
FCIQMC wiU also apply to DMQMC. Due to the relative 
youth of FCIQMC, the rate of theoretical and algorith- 
mic improvements is rapid [l6l . [2^ . 

The lattices studied with DMQMC so far are small 
by the standards of path-integral quantum Monte Carlo 
methods. However, the unique and defining feature of 
DMQMC is that it samples the full density matrix. Given 
the prominent role of the density matrix in the study of 
quantum information measures and the wide applicabil- 
ity of DMQMC, we hope that DMQMC wiU prove useful 
enough to justify further development. 
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with ay the Pauli spin matrix for the y-direction. This 
expression is only valid in the standard basis set. 

The value of the concurrence C ranges from zero to one 
and is monotonically related to the entanglement of for- 
mation [3^; the concurrence can therefore be regarded 
as a measure of entan glem ent. The entanglement of for- 
mation for two qubits [40| is given by 



£{C) = h 



l + ^/T^ 



where 



h{x) ^ -xlog2 X - (1 - x) log2(l - x). 



(A4) 



(A5) 



If the Hamiltonian and thus the reduced density matrix 
arc real (as is the case in the Heiscnberg model), the 
calculation of C reduces to the calculation of the moduli 
of the eigenvalues of 



R = PAia-y ® (Ty). 



(A6) 



Since i2 is a 4 x 4 matrix, it is trivial to compute the 
concurrence by direct diagonalization once pA is known. 

Using DMQMC, it is possible to sample the unnormal- 
ized reduced density matrix as described in Sec. IIIII and 
thus to estimate the concurrence using 



(C) = 



max(0, 71 - 72 - 73 - 74) 
Tr(pA) 



(A7) 



where pA is now the unnormalized reduced density ma- 
trix and {7i} are the eigenvalues of the unnormalized 
matrix R. 



Appendix B: Importance sampling 



The DMQMC evolution equation is 



(Bl) 



Appendix A: Concurrence 

For the special case of a two-qubit mixed state, the en- 
tanglement of formation can be calculated from a quan- 
tity known as the "concurrence" . Given a reduced den- 
sity matrix p^, where in this case A refers to a subsystem 
of two qubits, the concurrence is defined as 



C{pa) = max(0, 71 - 72 - 73 - 74), 



(Al) 



where 71 > 72 > 73 > 74 are the eigenvalues of the 
matrix 



R 



PaPa-sJPa 



(A2) 



and 



Pa = {Cy ® CFy)P*A{^y ® CTy), 



where pij is an element of p in the many-particle basis 
chosen for the simulation. Instead of p^j, we would like 
to sample the importance-sampled density matrix 



P13 



Ptj 



(B2) 



where E(i,j) is the excitation level of the pair («,j) and 
Wa is defined in Eq. (|2T|) . Following the standard proce- 
dure of importance sampling, we introduce a trial func- 
tion 



T 

P^J 



1 



Wt 



(B3) 



This matrix is symmetric, pfj = pj^, as E{i,j) = E{j,i). 
The importance-sampled density matrix then has com- 
(A3) ponents pij = pfjPij, where no summation is performed 



11 



over indices. Multiplying Eq. (|Bip by pfj yields the fol- 
lowing evolution equation for pij : 




The above differential equation is entirely analogous to 
Eq. (fT5|) . As such, the finite-difference version of Eq. (jB6p 



can be simulated in an almost identical manner to the 
standard DMQMC algorithm: the extra factors of 
simply act to alter the spawning probabilities. Consider 
the case where the excitation level of is 7 and ex- 
citation level of (fc, j) is 7 — 1. The probability that a 
spawning attempt from to is successful is al- 

tered by the factor 



Pkj ^7 

where Eq. ((2T|) has been used to simplify the expres- 
sion. Similarly, when spawning in the opposite direction, 
the probability is altered by the reciprocal of this factor. 
Spawning events that do not alter the excitation level are 
unaffected. 
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